% this function provide verification of partial P0: we start with BS model
% with const r and sigma, then the call price can be calculated as N(d1)
% N(d2) as in page 80 of math 133 book. so that the greeks can be computed
% exactly. see page 100 of math 133 book for analytical solutions.

function [P01,P02,delta1, delta2, gamma1, gamma2, pt3P01, pt3P02] ...
    = verify_partialP0

x = 4; r = 0.04; T = 1; t = 0; K=5; sigma = 0.3;

d1 = (log(x/K) + (r+1/2*sigma^2)*(T-t))/sigma/sqrt(T-t);
d2 = (log(x/K) + (r-1/2*sigma^2)*(T-t))/sigma/sqrt(T-t);

mu = r*(T-t)-sigma^2/2*(T-t); sigma = sqrt(sigma^2*(T-t));

P01 = x*normcdf(d1) - K*exp(-r*(T-t))*normcdf(d2);
P02 = exp(-r*T)*quadgk(@(y) max(y-K,0).*lognpdf(y/x,mu,sigma)/x,-Inf,Inf);

delta1 = normcdf(d1);
delta2 = exp(-r*T)*quadgk(@(y) max(y-K,0).*lognpdf(y/x,mu,sigma).*(log(y/x)-mu)/sigma^2/x^2, -Inf, Inf);

gamma1 = 1/sqrt(2*pi)*exp(-1/2*d1^2)/x/sigma/sqrt(T-t);
gamma2 = exp(-r*T)*quadgk(@(y) max(y-K,0).*lognpdf(y/x,mu,sigma).*((log(y/x)-mu).^2+sigma^2*(mu-log(y/x))-sigma^2)/sigma^4/x^3,-Inf, Inf);


pt3P01 = (1/sqrt(2*pi)*exp(-1/2*d1^2)*(-d1)-1/sqrt(2*pi)*exp(-1/2*d1^2)*sigma*sqrt(T-t))/...
    x^2/sigma^2/(T-t);
pt3P02 = exp(-r*T)*quadgk(@(y) max(y-K,0).*lognpdf(y/x,mu,sigma).*(-(sigma^2+mu-log(y/x)).*((mu-log(y/x)).^2+2*sigma^2*(mu-log(y/x))-3*sigma^2))/sigma^6/x^4,-Inf,Inf);
end